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*sO . Abstract The pair distribution function of monodisperse rigid spherocylinders is calculated by Shinomoto's 

method, which was originally proposed for hard spheres. The equation of state is derived by two different 
routes: Shinomoto's original route, in which a hard wall is introduced to estimate the pressure exerted on it, 
and the virial route. The pressure from Shinomoto's original route is valid only when the length-to-width ratio 
is less than or equal to 0.25 (i.e., when the spherocylinders are nearly spherical). The virial equation of state 
is shown to agree very well with the results of numerical simulations of spherocylinders with length-to-width 
ratio greater than or equal to 2. 
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1 Introduction 

■ Predicting the equation of state is a simple approach to probe the validity of approximation methods in the 

Q\ \ theory of liquids. For hard-sphere (HS) fluids, several methods have been proposed to derive the equation of 

state, such as the energy, virial, or compressibility equation of state in terms of an approximate pair distribution 
function; and the virial expansion of the equation of state [T|. Among the many derivations, the method 
proposed by Shinomoto is unique. 

Shinomoto derived the equation of state for HS fluids by calculating the density profile and pair dis- 
tribution function (2). The calculation is based on estimating the minimum work required to overcome the 
depletion force (3) between (i) an HS and a hard wall and (ii) between two HSs. Shinomoto's method is 
simpler and more concise than the numerous methods based on integral equations and is as accurate as the 
semiempirical Carnahan-Stirling equation of state 0J. Using Shinomoto's concept, the accuracy of the pair 
distribution function for an HS fluid turns out to be midway between that obtained by the Born-Green- Yvon 
theory and by the Percus-Yevick theory 01 . 

Shinomoto's method is not restricted to spherical molecules but is applicable to complex molecules as 
well. In this paper, I apply Shinomoto's concept to parallel hard-spherocylinder (SPC) fluids. By calculating 
the pair distribution function of the SPC fluid via Shinomoto's method, I obtain the fluid pressure by both 
Shinomoto's original route and the virial route [5 1. The results agree well with numerical simulation data via 
the former route when the length-to-width ratio is less than and equal to 0.25, and via the latter route for SPCs 
with length-to-width ratio greater than and equal to 2 up to °°. 
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2 Theory 

2.1 Shinomoto's Method for Hard-SPC Fluids 



In this section, I extend Shinomoto's method for a uniform monodisperse fluid consisting of parallel rigid 
SPCs. Although the discussion here is restricted to hard SPC fluids, it is also applicable, with slight modifi- 
cations, to general hard-body fluids. 

The length and width of the SPCs are represented by L and D, respectively, and the number density of 
the bulk SPC fluid is given by p. The packing fraction is denoted by 77 = pv, with the volume of an SPC 
being v = 7iD 3 /6 + nD 2 L/4. The symmetry axes of all SPCs in the fluid are aligned in a common direction. 
Since all SPCs are identical and parallel, the excluded volume of an SPC is an SPC with length 2L and width 
2D(KgUJi). 




Fig. 1 Schematic used to define excluded volume, (a) Each SPC (thick solid curve) has an excluded volume (inside dashed 
curve) into which the center of another SPC (thin solid curve) cannot enter, (b) Excluded volume (below the dashed line) of hard 
wall (shaded area). 



The SPC pair distribution function is 

*(r) = exp[-0*(r)], (1) 

where /3 is the inverse temperature and <£>(r) is the minimum work required to move a test particle (denoted 
by SPC t ) from infinity to point r under the condition in which another SPC (denoted by SPCo) is fixed at 
the origin of the coordinate system. The minimum work <J>(r) is nonzero because SPCt is subjected to a 
nonvanishing effective force F since the density pg(r) is now inhomogeneous because of the existence of 
SPCo- For example, when SPCo an d SPQ are close enough that their excluded volumes overlap, the pressure 
exerted on the surface of SPQ from surrounding SPCs vanishes inside the excluded volume of SPCo because 
of the absence of any SPC within this volume. However, outside this excluded volume, the pressure is nonzero. 
This mechanism to generate an effective force is similar to that required for the depletion force between 
particles suspended in the solutions of macromolecules (3J- 

Let S(t) be the surface of the excluded volume of SPQ at r, and let r* be a point on S(r). Since the kinetic 
pressure at r* is pg(r*)//3, the net force exerted on SPQ at r is 

F(r) = ^-/ g(r>(r*)dA(r*), (2) 

PV Js(r) 

where n(r*) is the inward unit vector normal to S at r* and dA(r*) is the infinitesimal area on S at r*. The 
integral is performed over S(r). The minimum work done against F is 



<&(r) = -^dr'-F(r'), (3) 
so the pair distribution function obeys the integral equation 



g(r) = exp 



- fV- f £(r*)n(r*)d4(r* 

V Joo Js(r') 



(4) 
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In principle, an approximate solution for this integral equation can be derived via an iteration method, as 
done in Ref. [4| for an HS fluid. In the present study, however, I follow the original paper by Shinomoto [2| 
and use the result of the first iteration. As an initial approximation to g(r), I adopt the low-density limit of 
the pair distribution function (i.e., exp[— j5u(r)], where u(r) denotes an intermolecular interaction between 
molecules separated by r). Since the low-density limit of g(r) is unity outside the excluded volume of SPCo 
and zero inside it, the integration over S on the right-hand side of Eq.© is zero if the SPC t excluded volume 
does not overlap with that of SPCo, an d nonzero if they do overlap. A simple calculation shows that the double 
integral on the right-hand side of Eq.© gives the overlap volume of the excluded volumes of SPC t and SPCo- 
Thus, 

g(r)=exp[Tjy ov:s (r)/v], (5) 

where V v:s(r) is the overlap volume of the excluded volumes of two SPCs separated by r. Equation (|5]l is 
alternatively expressed as exp[r/ / dr,/o;/(t] in terms of Mayer's /-function /q,- between SPCo an d an SPC 
labeled i, and / t between i and SPC t . The expansion of Eq.© in powers of 77 gives an /-bond expansion of 
g(r) which is exact up to the first order in 77. 

To calculate the equation of state as the pressure exerted on a wall, I introduce a hard wall at z = 0. Let the 
angle between the SPC symmetry axes and the z axis be 9, as shown in FigJTJ). An SPC center cannot enter the 
excluded volume of the wall; that is, the SPC center cannot enter the region z < z c , where z c = (D+Lcos 9)/2 
is the z coordinate of the point of contact. To derive the first approximation of the density profile, consider the 
minimum work to move SPC t toward the hard wall (a derivation similar to that of Eq.©), which gives 

pi(r,e) = pi(z,e) = ^exp[77V ov:w (z,0)/v], (6) 

where V ov - w (z, 6) is the overlap volume of the excluded volumes of the wall and SPC at r = (x,y,z). 

The pair distribution function I© and the density approximation I© improve the approximation of the pres- 
sure on the test particle and thus the minimum work. Since the density at point r* on S(r) is p\ (r*, 9)g(r* — r), 
the pressure at r* takes the form 

px (r*,r,6) = pi (r*, 9)g(r* -r)/j8 = exp{T] [V ov:s (r* - r) + V ov:w (z*, 9)}/v}. (7) 

The minimum work to move a test SPC from infinity to the point of contact with the wall is 

<P c (d) = - T dz' [ Pl (r*y,e)n z (r*)dA(r*), (8) 

where S'(r) denotes the subset of 5(r) outside the excluded volume of the wall. Again, the minimum work © 
defines the contact density at the wall, which is pexp[— /3<£> c ]. Applying the exact pressure sum rule [6] leads 
to the conclusion that the dimensionless bulk pressure P* = [iPv is equivalent to the contact packing fraction 
at the wall, so that the resulting dimensionless pressure is 

P s *(e) = T]e X p[-/3<& c (e)]. (9) 

The subscript S indicates that, when using Eq.dTOll to estimate /3<P C (9), the pressure is the result of a straight- 
forward extension of Shinomoto's method to parallel SPCs. 

By expanding pi (r* , r', 9) in powers of the packing fraction T] up to the second order and substituting the 
result into Eq.®, we obtain the following power expansion of P<P C (9)'- 

j3* c (0) = -4T]-T] 2 4 TV/ [Wr*-O+V ov:w (^0)K(r*)dA(r*/ (10) 

V L J 00 Js'(r) 

The coefficient of the first term is independent of L and 9 since it is given by — V v:w(zc, 0)/v, where 
Vov:w(zc, ^) i s na lf °f an excluded volume for any 9. Note that the coefficient of the first term in Eq.lfTOil 
gives the exact second virial coefficient (7j if Eq.© is expanded in powers of 77. Since Eq. dlOl) contains Vov:w 
the pair distribution function exp[j34> c ] cannot be compared to the usual /-bond expansion of g(r), as done 
below Eq.©. 

It is possible to obtain V v:s(r) and V ov - v/ (z,9) as combinations of the overlap volumes of hemispheres, 
cylinders, and a half-infinite region. Here I omit explicit expressions for these functions since the calculation 
is elementary but the resulting expressions are considerably tedious. 
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The second term in the surface integral in Eq. dlOb can be reduced to integration over the single variable 
Z*, since V v:w(z* , 0) is constant when z* is fixed, so the integration over the other variable can be performed 
analytically. Similarly, integration of the first term in the surface integral can be reduced to a single variable 
integration along the SPC axis. The remaining double integral in Eq. dlOl) was calculated numerically. 



2.2 Virial Equation of State 

In this section, I derive the equation of state through the virial route [5] (i.e., not strictly following Shinomoto's 
method as in the previous subsection) by using the pair distribution function (Eq.©) without introducing a 
hard wall. The virial equation of state is 



(11) 



where u(r) denotes an intermolecular interaction between molecules separated by r. In our case, Vw(r) van- 
ishes unless r corresponds to a contact configuration, so only the contact value of g(r) contributes to the virial 
pressure. 

In the molecular frame, in which the z axis lies along the SPC symmetry axis, the contact value of g(r) = 
g(x,y,z) is a function of the single variable z. This may be seen by noting that, although g(r) is a function 
of z and x 2 + y 2 , the sum xr + y 2 is defined by z provided that r corresponds to a contact configuration. The 
contact value of g(r) in Eq.© is denoted by g c (z) ■ The dimensionless virial pressure thus reduces to 



27TT] 2 D 3 
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K = n + 



- rL A Z nL+\( {Z-L)L\ 



8e(z) 



(12) 



3 Results 

With isobaric Monte Carlo (MC) simulations [8], I obtained "exact" results for g c (z) and the equation of 
state, which are to be compared below with Shinomoto and virial analytical results. Note that the following 
MC results are not used in the evaluations of P| and P* in the previous section. In the MC simulations, I 
used a length-to-width ratio L/D = 5 and N = 4096 SPCs. The SPCs are perfectly aligned along the z axis 
in a simulation box, and the periodic boundary conditions are imposed in all directions. I performed 10 6 MC 
steps and recorded g c (z) and the packing fraction every 10 3 MC steps to calculate their average values. For 
isobaric MC simulations, I also checked the consistency of the given pressure by comparing it to the pressure 
calculated with the right-hand side of Eq. dl2l) with the pair distribution function calculated via an isochoric 
MC simulation (the pressures agree within the error bars). For L/D = 5 parallel SPCs, MC simulations [9| 
give a pressure and the packing fraction at the nematic-smectic phase transition off* = 2.1 1 and 77 = 0.406, 
respectively, for a system of N = 270 SPCs. The results of molecular dynamics simulations [10] are the same 
for a system of N — 1080 SPCs. Since the discussion in the previous section is valid only for the nematic 
phase, I performed MC simulations at P* = 0.5, 1.0, 1.5, and 2.0; under these conditions, no distinct layer 
structure was found. 

Figure [2] shows the contact values of g c (z) given by Eq.© and those obtained from the MC simulations. At 
low density (77 = 0.2076), Eq.© agrees well with the contact value of the pair distribution function obtained 
from the isochoric MC simulation. However, the results diverge at high density. For example, for 2.8 < z/D < 
5> gc(z) for T) = 0.3929 is lower than g c (z) for r\ = 0.2076; such an inverse dependency on 77 cannot be 
expected on the basis of Eq.((5]l since g c (z) in Eq.© increases monotonically as a function of V\ for a given 
Z. A similar failure of Shinomoto's pair distribution function at high density has also been reported for rigid 
parallel cylinders ifTTTl . 

The agreement between the equation of state and the simulation data is all the more striking given the 
poor accuracy of our estimate for the pair distribution function at high density. Figure [3] shows the equation 
of state for SPCs with L/D = 5. The virial pressure P* agrees very well with the results of MC simulations 
throughout the nematic phase. The empirically modified result of the scaled particle theory for parallel hard 
SPCs [12], P S * PT = (2?7 2 + ?])/(l - T]) 2 , is also shown in Fig[3j P* agrees better with MC simulation results 

than PXryy. 
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Fig. 2 Contact value of pair distribution function versus z/D (the normalized SPC separation in the symmetry axis direction) for 
SPCs with L/D = 5. The open and solid circles correspond to isochoric MC simulations for packing fractions T] = 0.2076 (for 
P* = 0.5) and T] = 0.3929 (for P* = 2.0), respectively. The solid and dashed curves indicate the respective contact values from 
Eq.l(5} for the same packing fractions. 




Fig. 3 Equations of state giving the pressure from Shinomoto's original route and the virial pressure as functions of packing 
fraction r\ for the nematic phase of parallel SPCs with L/D = 5. The solid and dashed curves indicate P| and P*, respectively. 
The crosses and circles indicate, respectively, data from isobaric and isochoric MC simulations. Their error bars are smaller than 
the symbol size. The dotted curve is the empirically modified result of the scaled particle theory in Ref, fl2l . 

The numerical evaluation of Eq. dlOl) indicates that <P C (9), and thus the pressure P$(9), are independent 
of 9. Although this result does not constitute a rigorous proof of the independence of these quantities with 
respect to 9, it does mean that the 9 -dependence of these quantities can be safely neglected. The pressure 
Pg , which agrees well with the exact result for HSs [2|, deviates from the results of MC simulations at high 
density for SPCs with L/D = 5. 

The dependence of P* on L/D for parallel SPCs is quite weak. The virial pressure with L/D = (i.e., for 
HSs) and L/D = °° is shown in Fig. |4j For L/D = 0, the exact expression for the virial pressure becomes 



p; = T ]+477 2 e 5r '/ 2 . 



(13) 
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The difference in virial pressures between SPCs with L/D — and °° is a few percent at T] = 0.444, which is 




0.1 0.2 0.3 0.4 0.5 
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Fig. 4 Equations of state giving the pressure from Shinomoto's original route and the virial pressure as functions of packing 
fraction r\ for the nematic phase of parallel SPCs with different length-to-width ratios. The solid and dashed curves give the 
virial pressure P* for SPCs with L/D = °° and (i.e., hard spheres), respectively. The dotted curve represents the pressure P$ 
for SPCs with L/D = 0.25. The squares, circles, triangles, and diamonds show the simulation results for SPCs with L/D = °°, 2, 
1, and 0.25, respectively. These simulation data are taken from Ref. 1 10| for L/D = <x and from Ref. (9) for the others. 

the packing fraction of the nematic-smectic phase transition for SPCs with L/D — °° ifTOl . All virial pressures 
for SPCs with finite L/D are between these two curves (see solid and dashed curves in Figj4j. 

The simulation data in Fig.[4]are taken from the molecular dynamics simulations for N = 1080 SPCs with 
L/D = oo (TO), and MC simulations for N = 90 SPCs with L/D = 2, 1 and 0.25 @. The weak dependence of 
the pressure on L/D is also found in numerical simulations if L/D > 2; thus, P* agrees well with these data 
from L/D = 2 to °°. For SPCs with L/D < 2, however, the pressure at a given 77 increases as L decreases, and 
the agreement between the virial pressure P* and the numerical results worsens. As expected from the success 
of Shinomoto's method for HSs, the pressure PS gives better results for SPCs with small L/D; but its validity 
is limited to < L/D < 0.25. 

4 Discussion 

The failure of g c (z) at high densities indicates that the agreement between virial pressures for SPCs from 
L/D = 2 to L/D = 00 and those of MC simulations is due to a coincidental cancellation of errors. Further- 
more, the virial method for calculating pressure is not valid for SPCs with small L/D, so we are obliged to 
use Shinomoto's original method; however, there is no theoretical justification to use two different methods 
depending on L/D. Yet, despite these flaws in the theory, it is interesting that P* agrees well with the numeri- 
cal data over a broad range of L/D, and that P| for L/D = and P* for L/D = 00 seem to give the upper and 
lower bounds for pressures for all L/D. 

In principle, Shinomoto's method can be generalized to a system of SPCs with arbitrary direction and 
even to arbitrary hard-core systems. However, it is difficult to implement this method for such systems. For 
example, in many cases, the excluded volume has a complex shape, so calculating the overlap volume of 
V exc (I2,I2") and V exc (Q" ,Q') is almost impossible (where Q is the orientation of a nonspherical molecule 
and V exc (Q , Q') is the excluded volume of two molecules with orientations Q and Q.'). Numerical estimation 
of the overlap volume might be possible by MC integration, as done in calculating the virial coefficients 
of nonspherical molecules. If such a numerical estimation is performed, Shinomoto's method, which is free 
from divergence, would have an advantage to the virial expansion, which has a considerably small radius of 
convergence for strongly anisotropic molecules flBl . 
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There is no explanation for the insensitivity of the equation of state to the SPC length in 2 < L/D < °°. 
If the system is a fluid of parallel hard ellipsoids, the dimensionless pressure is independent of the molecular 
length since the system can be mapped to an HS fluid with a scale transformation. However, the SPC fluid 
cannot be mapped onto a single system and the above discussion for ellipsoids is not applicable. 



5 Conclusion 

The pair distribution function for an SPC fluid derived by Shinomoto's method agrees well with the result 
of MC simulations at low density but the agreement worsens at high densities. From this Shinomoto's pair 
distribution function, the equation of state is derived by two different routes: Shinomoto's original route and 
the virial route. Despite the failure of Shinomoto's pair distribution function, the virial pressure P* agrees very 
well with the numerical data for 2 < L/D < °° throughout the nematic phase. The pressure from Shinomoto's 
original route, P|, agrees well with the numerical data for < L/D < 0.25. P| for L/D = and P* for 
L/D — > °° seem to give the upper and lower bounds for pressures for all L/D. 
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